library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.1
library(psych) #for pairs.panels, but could use other packages, e.g. GGalley
## Warning: package 'psych' was built under R version 4.1.1
library(lavaan)
## Warning: package 'lavaan' was built under R version 4.1.1
library(semPlot)
## Warning: package 'semPlot' was built under R version 4.1.1
library(DiagrammeR)
## Warning: package 'DiagrammeR' was built under R version 4.1.1

Import data

combined=read.csv("data/annual_averages/annual_data_compiled.csv")
cnames=read.csv("analysis/column_names.csv", stringsAsFactors = F)
dsub=filter(combined, Year>=1975)
focaldata=dsub[,cnames$Datacolumn]
fvars=cnames$Shortname
colnames(focaldata)=fvars

source("analysis/myLavaanPlot.r")

Data prep

Log transform, scale

#log transform
logvars=fvars[cnames$Log=="yes"]
logtrans=function(x) {
  x2=x[which(!is.na(x))]
  if(any(x2==0)) {log(x+min(x2[which(x2>0)],na.rm=T))}
  else {log(x)}
}
focaldatalog = focaldata %>% 
  mutate_at(logvars,logtrans)

#scale data
fd=focaldatalog
fd=fd %>% 
  mutate_at(2:length(fvars),list("1"=lag)) %>% #lag 1
  mutate_at(2:length(fvars),list("fd"=function(x) c(NA,diff(x)))) %>% #first difference
  mutate_at(2:length(fvars),list("dtr"=function(x) { #detrend
    x2=x
    x2[x2==0]=NA
    res=residuals(lm(x2~fd$year))
    out=x
    out[which(!is.na(x2))]=res
    return(out)
  })) %>%
  mutate_at(-1,scale)

Time series plots

Original units

Log scaled

First difference

Detrended

Bivariate plots

SEM model

Best model so far

model4='zoop=~hcope+clad+mysid
        fish=~estfish+estfish_bsmt+estfish_bsot
        zoop~sside+chla
        chla~clams+flow
        fish~zoop+flow+sside
        hcope~~mysid
        clad~~mysid
        hcope~~clad'
modfit4=sem(model4, data=fd)
summary(modfit4, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 39 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        23
##                                                       
##                                                   Used       Total
##   Number of observations                            40          46
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                34.739
##   Degrees of freedom                                26
##   P-value (Chi-square)                           0.117
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   zoop =~                                                               
##     hcope             1.000                               0.607    0.615
##     clad              1.341    0.334    4.013    0.000    0.814    0.913
##     mysid             1.421    0.318    4.461    0.000    0.862    0.921
##   fish =~                                                               
##     estfish           1.000                               0.758    0.766
##     estfish_bsmt      1.153    0.197    5.852    0.000    0.874    0.899
##     estfish_bsot      0.968    0.193    5.011    0.000    0.734    0.772
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   zoop ~                                                                
##     sside            -0.513    0.155   -3.298    0.001   -0.844   -0.515
##     chla              0.409    0.123    3.334    0.001    0.674    0.531
##   chla ~                                                                
##     clams            -0.415    0.112   -3.689    0.000   -0.415   -0.495
##     flow             -0.342    0.106   -3.215    0.001   -0.342   -0.431
##   fish ~                                                                
##     zoop              0.558    0.315    1.770    0.077    0.447    0.447
##     flow              0.303    0.100    3.041    0.002    0.400    0.397
##     sside            -0.391    0.240   -1.628    0.104   -0.516   -0.315
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .hcope ~~                                                              
##    .mysid             0.103    0.145    0.708    0.479    0.103    0.362
##  .clad ~~                                                               
##    .mysid            -0.097    0.154   -0.628    0.530   -0.097   -0.729
##  .hcope ~~                                                              
##    .clad              0.019    0.134    0.141    0.888    0.019    0.067
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .hcope             0.605    0.183    3.298    0.001    0.605    0.621
##    .clad              0.132    0.157    0.843    0.399    0.132    0.167
##    .mysid             0.134    0.174    0.766    0.444    0.134    0.152
##    .estfish           0.405    0.108    3.759    0.000    0.405    0.413
##    .estfish_bsmt      0.181    0.077    2.345    0.019    0.181    0.192
##    .estfish_bsot      0.365    0.098    3.729    0.000    0.365    0.404
##    .chla              0.418    0.093    4.472    0.000    0.418    0.675
##    .zoop              0.123    0.103    1.197    0.231    0.334    0.334
##    .fish              0.171    0.076    2.235    0.025    0.297    0.297
## 
## R-Square:
##                    Estimate
##     hcope             0.379
##     clad              0.833
##     mysid             0.848
##     estfish           0.587
##     estfish_bsmt      0.808
##     estfish_bsot      0.596
##     chla              0.325
##     zoop              0.666
##     fish              0.703
semPaths(modfit4, "std", edge.label.cex = 1, residuals = F)

semPaths(modfit4, "par", edge.label.cex = 1, residuals = F)

# lavaan::parameterestimates(modfit4)

labels1 <- createLabels(modfit4, cnames)

## Diagram without model covariances:
myLavaanPlot(model=modfit4, labels=labels1, 
                         node_options=list(shape="box", fontname="Helvetica"), 
                         coefs=TRUE, stand=TRUE, covs=FALSE, sig=0.05, 
                         width=c("regress","latent"),
                         color=c("regress","latent"))
## Diagram with model covariances:
myLavaanPlot(model=modfit4, labels=labels1, 
                         node_options=list(shape="box", fontname="Helvetica"), 
                         coefs=TRUE, stand=TRUE, covs=TRUE, sig=0.05, 
                         width=c("regress","latent","covs"),
                         color=c("regress","latent","covs"))

Same model with detrended data

model4b='zoop=~hcope_dtr+clad_dtr+mysid_dtr
        fish=~estfish_dtr+estfish_bsmt_dtr+estfish_bsot_dtr
        zoop~sside_dtr+chla_dtr
        chla_dtr~clams_dtr+flow_dtr
        fish~zoop+flow_dtr+sside_dtr'
modfit4b=sem(model4b, data=fd)
summary(modfit4b, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 29 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        20
##                                                       
##                                                   Used       Total
##   Number of observations                            40          46
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                34.591
##   Degrees of freedom                                29
##   P-value (Chi-square)                           0.218
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   zoop =~                                                               
##     hcope_dtr         1.000                               0.743    0.754
##     clad_dtr          0.867    0.245    3.543    0.000    0.644    0.641
##     mysid_dtr         1.072    0.269    3.990    0.000    0.796    0.795
##   fish =~                                                               
##     estfish_dtr       1.000                               0.441    0.485
##     estfsh_bsmt_dt    2.169    0.769    2.818    0.005    0.956    0.971
##     estfsh_bst_dtr    1.296    0.482    2.688    0.007    0.571    0.580
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   zoop ~                                                                
##     sside_dtr        -0.303    0.182   -1.670    0.095   -0.408   -0.272
##     chla_dtr          0.390    0.138    2.832    0.005    0.525    0.492
##   chla_dtr ~                                                            
##     clams_dtr         0.056    0.128    0.435    0.664    0.056    0.063
##     flow_dtr         -0.411    0.137   -3.007    0.003   -0.411   -0.435
##   fish ~                                                                
##     zoop              0.104    0.100    1.036    0.300    0.175    0.175
##     flow_dtr          0.248    0.106    2.348    0.019    0.563    0.558
##     sside_dtr        -0.093    0.099   -0.946    0.344   -0.212   -0.141
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .hcope_dtr         0.419    0.145    2.892    0.004    0.419    0.431
##    .clad_dtr          0.595    0.161    3.685    0.000    0.595    0.590
##    .mysid_dtr         0.369    0.149    2.467    0.014    0.369    0.368
##    .estfish_dtr       0.633    0.148    4.267    0.000    0.633    0.765
##    .estfsh_bsmt_dt    0.055    0.191    0.291    0.771    0.055    0.057
##    .estfsh_bst_dtr    0.644    0.160    4.017    0.000    0.644    0.664
##    .chla_dtr          0.697    0.156    4.472    0.000    0.697    0.795
##    .zoop              0.395    0.168    2.352    0.019    0.716    0.716
##    .fish              0.121    0.077    1.573    0.116    0.620    0.620
## 
## R-Square:
##                    Estimate
##     hcope_dtr         0.569
##     clad_dtr          0.410
##     mysid_dtr         0.632
##     estfish_dtr       0.235
##     estfsh_bsmt_dt    0.943
##     estfsh_bst_dtr    0.336
##     chla_dtr          0.205
##     zoop              0.284
##     fish              0.380
semPaths(modfit4b, "std", edge.label.cex = 1, residuals = F)

semPaths(modfit4b, "par", edge.label.cex = 1, residuals = F)

Old models (not run)

model1='clams~year+flow
        sside~year+flow
        chla~flow+clams #+year?
        hcope~chla+pcope #+year?
        clad~chla+pcope
        #pcope~year+clams
        mysid~hcope+pcope
        estfish~hcope+pcope+mysid+sside'
modfit1=sem(model1, data=fd)
summary(modfit1, standardized=T, rsq=T)
semPaths(modfit1, "std", edge.label.cex = 1, residuals = F)

residuals(modfit1)

model1='clams~year+flow
        sside~year+flow
        foodsupply=~chla+hcope+clad+mysid
        foodsupply~clams+flow+sside
        fish=~estfish+estfish_bsmt
        fish~foodsupply'
modfit1=sem(model1, data=fd)
summary(modfit1, standardized=T, rsq=T)
semPaths(modfit1, "std", edge.label.cex = 1, residuals = F)

model1='foodsupply=~chla+hcope+clad+mysid
        fish=~estfish+estfish_bsmt
        foodsupply~clams+flow+sside
        fish~foodsupply'
modfit1=sem(model1, data=fd)
summary(modfit1, standardized=T, rsq=T)
semPaths(modfit1, "std", edge.label.cex = 1, residuals = F)

model2='foodsupply=~chla+hcope+clad+mysid
        fish=~estfish+estfish_bsmt
        foodsupply~clams+flow+sside
        fish~foodsupply+flow+sside'
modfit2=sem(model2, data=fd)
summary(modfit2, standardized=T, rsq=T)
semPaths(modfit2, "std", edge.label.cex = 1, residuals = F)

residuals(modfit2)


model3='zoop=~hcope+clad+mysid
        fish=~estfish+estfish_bsmt
        zoop~clams+flow+sside+chla
        chla~clams #+tphos+nitrate+ammonia
        fish~zoop+flow+sside'
modfit3=sem(model3, data=fd)
summary(modfit3, standardized=T, rsq=T)
semPaths(modfit3, "std", edge.label.cex = 1, residuals = F)

residuals(modfit3, type="cor")
modificationindices(modfit3)

model4='zoop=~hcope+clad+mysid
        fish=~estfish+estfish_bsmt+estfish_bsot
        zoop~sside+chla
        chla~clams+flow
        fish~zoop+flow+sside
        hcope~~mysid
        clad~~mysid
        hcope~~clad'
modfit4=sem(model4, data=fd)
summary(modfit4, standardized=T, rsq=T)
semPaths(modfit4, "std", edge.label.cex = 1, residuals = F)
semPaths(modfit4, "par", edge.label.cex = 1, residuals = F)

residuals(modfit4, type="cor")
modificationindices(modfit4)
#using estfish instead of smelt
model1f='chla~year+flow+secchi
        hcope~year+flow+temp+chla+pcope+secchi
        pcope~year
        mysid~year+hcope+pcope
        estfish~year+secchi+hcope+pcope+mysid'
modfit1f=sem(model1f, data=fd)
#summary(modfit1f, standardized=T, rsq=T)
semPaths(modfit1f, "std", edge.label.cex = 1, residuals = F)

#using marfish instead of smelt
model1g='chla~year+flow+secchi
        hcope~year+flow+temp+chla+pcope+secchi
        pcope~year
        mysid~year+hcope+pcope
        marfish~year+secchi+hcope+pcope+mysid'
modfit1g=sem(model1g, data=fd)
#summary(modfit1g, standardized=T, rsq=T)
semPaths(modfit1g, "std", edge.label.cex = 1, residuals = F)

Detrended (anomalies) without year effect

#smelt
model2a='chla_dtr~flow_dtr+temp_dtr+secchi_dtr
         hzoop_dtr~flow_dtr+temp_dtr+secchi_dtr+chla_dtr+pcope_dtr
         pcope_dtr~flow_dtr+temp_dtr+secchi_dtr
         mysid_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr
         smelt_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr+mysid_dtr'
modfit2a=sem(model2a, data=fd)
#summary(modfit2a, standardized=T, rsq=T)
semPaths(modfit2a, "std", edge.label.cex = 1, residuals = F)

#estuarine fishes
model2b='chla_dtr~flow_dtr+temp_dtr+secchi_dtr
         hzoop_dtr~flow_dtr+temp_dtr+secchi_dtr+chla_dtr+pcope_dtr
         pcope_dtr~flow_dtr+temp_dtr+secchi_dtr
         mysid_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr
         estfish_dtr~flow_dtr+temp_dtr+secchi_dtr+hzoop_dtr+pcope_dtr+mysid_dtr'
modfit2b=sem(model2b, data=fd)
#summary(modfit2b, standardized=T, rsq=T)
semPaths(modfit2b, "std", edge.label.cex = 1, residuals = F)